# Making univariate toy data
library(dplyr,
quietly = TRUE)
mz_data <- MASS::mvrnorm(
n = 100,
mu = c(0, 0),
Sigma = matrix(c(1.0, 0.8,
0.8, 1.0),
nrow = 2)) %>%
as.data.frame() %>%
setNames(c('X1', 'X2'))
dz_data <- MASS::mvrnorm(
n = 100,
mu = c(0, 0),
Sigma = matrix(c(1.0, 0.5,
0.5, 1.0),
nrow = 2)) %>%
as.data.frame() %>%
setNames(c('X1', 'X2'))
data <- rbind(data.frame(zyg = 'MZ', mz_data),
data.frame(zyg = 'DZ', dz_data)) %>%
mutate(
zyg = factor(zyg, levels = c('MZ', 'DZ'))
)
head(data)
## zyg X1 X2
## 1 MZ 0.55107062 -0.54415443
## 2 MZ -1.77073118 -1.81837605
## 3 MZ 1.45332458 1.40457435
## 4 MZ -0.28017613 0.17413459
## 5 MZ -0.55666668 -1.33562604
## 6 MZ -0.03354105 0.01393325# Installing the packages
# install.packages('devtools')
# library(devtools)
# install_github('ivanvoronin/mlth.data.frame')
# install_github('ivanvoronin/TwinAnalysis')
# Loading the packages
library(OpenMx)
library(TwinAnalysis)
# Descriptive statistics
get_univ_descriptives(data, zyg = 'zyg', vars = 'X')
## $X
## Twin1-------------------- Twin2------------------- Cor--------------------------
## n M SD n M SD r lowerCI upperCI
## MZ 100 0.10347156 1.0247719 100 0.02369611 1.015420 0.8394995 0.7700587 0.8892814
## DZ 100 -0.01673161 0.9709349 100 -0.14933507 1.038595 0.5738636 0.4253948 0.6922529# Univariate ACE model
model <- univariate_ace(data, zyg = 'zyg', ph = 'X')
model <- mxRun(model)
summary(model)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 -0.0154077 0.06578862
## 2 ACE.a[1,1] a X X 0.7347096 0.08905736
## 3 ACE.c[1,1] c X X 0.5603042 0.12403442
## 4 ACE.e[1,1] e X X 0.4102589 0.02890131
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 4 396
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 984.723
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/400
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 192.723 992.723
## BIC: -1113.411 1005.916
## | Sample-Size Adjusted
## AIC: 992.9281
## BIC: 993.2438
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:50:48
## Wall clock time: 0.07298517 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)# Constrained models
# This is a list of MxModels
nested_models <- def_nested_models(
model,
AE = mxConstraint(c[1, 1] == 0),
CE = mxConstraint(a[1, 1] == 0),
E = list(mxConstraint(c[1, 1] == 0),
mxConstraint(a[1, 1] == 0)),
run = TRUE
)
fit_stats(model, nested_models = nested_models)
## Warning in computeFitStatistics(likelihood, DoF, chi,
## chiDoF, retval[["numObs"]], : Your model may be mis-
## specified (and fit worse than an independence model),
## or you may be using the wrong independence model, see ?
## mxRefModels
## model base ep minus2LL df AIC BIC
## 1 ACE Saturated 4 984.7230 396 192.7230 -1113.4107
## 2 AE ACE 4 989.1442 397 195.1442 -1114.2878
## 3 CE ACE 4 1007.3842 397 213.3842 -1096.0478
## 4 E ACE 4 1144.5469 398 348.5469 -964.1834
## CFI TLI RMSEA diffLL diffdf
## 1 1.0005878 1.0001959 0.00000000 5.905993 6
## 2 0.9791975 0.9940564 0.04875019 4.421220 1
## 3 0.8651567 0.9614733 0.12411747 22.661211 1
## 4 0.0138350 0.7534587 0.31397644 159.823929 2
## p
## 1 4.338028e-01
## 2 3.549466e-02
## 3 1.932294e-06
## 4 1.970946e-35# Making bivariate toy data
mz_data2 <- MASS::mvrnorm(
n = 100,
mu = c(0, 0, 0, 0),
Sigma = matrix(
c(1.0, 0.7, 0.6, 0.5,
0.7, 1.0, 0.5, 0.8,
0.6, 0.5, 1.0, 0.7,
0.5, 0.8, 0.7, 1.0),
nrow = 4)) %>%
as.data.frame() %>%
setNames(c('X1', 'Y1', 'X2', 'Y2'))
dz_data2 <- MASS::mvrnorm(
n = 100,
mu = c(0, 0, 0, 0),
Sigma = matrix(
c(1.0, 0.7, 0.4, 0.2,
0.7, 1.0, 0.2, 0.4,
0.4, 0.2, 1.0, 0.7,
0.2, 0.4, 0.7, 1.0),
nrow = 4)) %>%
as.data.frame() %>%
setNames(c('X1', 'Y1', 'X2', 'Y2'))
data2 <- rbind(data.frame(zyg = 'MZ', mz_data2),
data.frame(zyg = 'DZ', dz_data2)) %>%
mutate(
zyg = factor(zyg, levels = c('MZ', 'DZ'))
)
head(data2)
## zyg X1 Y1 X2 Y2
## 1 MZ 2.9737961 1.7046050 2.6551091 0.93346789
## 2 MZ -0.2515991 -0.1794637 -0.8486091 -0.65641665
## 3 MZ -1.2586533 -0.5918007 -1.1710607 -0.08288749
## 4 MZ 0.4271026 0.3153837 0.6676022 0.78442769
## 5 MZ 0.2249155 0.5210256 -0.4747569 -1.04644136
## 6 MZ 0.3755235 0.4408710 -0.4999070 -0.53263830# Descriptive statistics
get_univ_descriptives(data2,
zyg = 'zyg',
vars = c('X', 'Y'))
## $X
## Twin1------------------ Twin2------------------- Cor--------------------------
## n M SD n M SD r lowerCI upperCI
## MZ 100 0.06760087 1.146399 100 0.02576845 0.9835557 0.6359860 0.5023170 0.7399680
## DZ 100 0.01026312 1.000688 100 0.01358189 1.1121019 0.4952969 0.3310888 0.6303876
##
## $Y
## Twin1-------------------- Twin2-------------------- Cor--------------------------
## n M SD n M SD r lowerCI upperCI
## MZ 100 -0.01586990 0.9520076 100 0.02848424 0.9729878 0.7641308 0.6679673 0.8351952
## DZ 100 -0.01175412 1.0119057 100 -0.02396046 1.0687223 0.5158921 0.3554967 0.6467724# Cross-twin cross-trait correlations
get_cross_trait_cors(data2,
zyg = 'zyg',
vars = c('X', 'Y'))
## zyg: MZ
## X1 Y1 X2 Y2
## X1 1.0000000 0.7778976 0.6359860 0.4868070
## Y1 0.7778976 1.0000000 0.5482844 0.7641308
## X2 0.6359860 0.5482844 1.0000000 0.7080632
## Y2 0.4868070 0.7641308 0.7080632 1.0000000
## --------------------------------------------
## zyg: DZ
## X1 Y1 X2 Y2
## X1 1.0000000 0.7242958 0.4952969 0.3904858
## Y1 0.7242958 1.0000000 0.2686450 0.5158921
## X2 0.4952969 0.2686450 1.0000000 0.7006547
## Y2 0.3904858 0.5158921 0.7006547 1.0000000# Bivariate ACE model
model2 <- multivariate_ace(
data = data2,
zyg = 'zyg',
ph = c('X', 'Y')
)
model2 <- mxRun(model2)
summary(model2)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 0.028582736 0.06602367
## 2 ACE.mean[1,2] mean 1 2 -0.006697648 0.06398131
## 3 ACE.a[1,1] a X X 0.229633789 0.19036365
## 4 ACE.a[2,1] a Y X 0.487159593 0.14552821
## 5 ACE.a[2,2] a Y Y 0.761979977 0.09835996
## 6 ACE.c[1,1] c X X 0.525658510 0.07115624
## 7 ACE.c[2,1] c Y X 0.367411377 0.20938628
## 8 ACE.c[2,2] c Y Y 0.456055458 0.15872435
## 9 ACE.e[1,1] e X X 0.417842833 0.02948707
## 10 ACE.e[2,1] e Y X 0.493955852 0.05375753
## 11 ACE.e[2,2] e Y Y 0.467266109 0.03298510
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 11 789
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 1793.504
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/800
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 215.5039 1815.504
## BIC: -2386.8685 1851.785
## | Sample-Size Adjusted
## AIC: 1816.908
## BIC: 1816.936
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:50:49
## Wall clock time: 0.077667 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)# Output tables from the bivariate ACE
get_output_tables(model2)
## $Variance_components
## A(%) C(%) E(%)
## X 0.2589905 0.3672559 0.3737536
## Y 0.5766131 0.2065536 0.2168333
##
## $Covariation_components
## Total A C E A(%)
## X <-> X 1.1199492 0.2900561 0.4113080 0.4185850 0.2589905
## X <-> Y 0.7695746 0.3712059 0.1675600 0.2308088 0.4823520
## Y <-> Y 1.0069377 0.5806135 0.2079866 0.2183376 0.5766131
## C(%) E(%)
## X <-> X 0.3672559 0.3737536
## X <-> Y 0.2177306 0.2999174
## Y <-> Y 0.2065536 0.2168333
##
## $All_correlations
## r_A r_C r_E Total
## X <-> X 1.0000000 1.0000000 1.0000000 1.0000000
## X <-> Y 0.9045451 0.5728871 0.7634776 0.7246867
## Y <-> Y 1.0000000 1.0000000 1.0000000 1.0000000A model (mxModel):
# Twin models in OpenMx
library(OpenMx)
# Univariate ACE model
model <- mxModel(
name = 'ACE',
# Means
mxMatrix(type = 'Full',
nrow = 1, ncol = 1,
free = TRUE,
name = 'mean'),
mxAlgebra(cbind(mean, mean),
name = 'expMeans'),
# Variance components
mxMatrix(type = 'Lower',
nrow = 1, ncol = 1,
free = TRUE,
name = 'a'),
mxMatrix(type = 'Lower',
nrow = 1, ncol = 1,
free = TRUE,
name = 'c'),
mxMatrix(type = 'Lower',
nrow = 1, ncol = 1,
free = TRUE,
name = 'e'),
mxAlgebra(t(a) %*% a, name = 'A'),
mxAlgebra(t(c) %*% c, name = 'C'),
mxAlgebra(t(e) %*% e, name = 'E'),
# Expected covariance
mxAlgebra(rbind(cbind(A + C + E, A + C ),
cbind(A + C, A + C + E)),
name = 'expCovMZ'),
mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
cbind(.5%x%A + C, A + C + E)),
name = 'expCovDZ'),
# MZ submodel
mxModel(name = 'MZ',
mxData(observed = mz_data,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovMZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'X2')),
mxFitFunctionML()),
# DZ submodel
mxModel(name = 'DZ',
mxData(observed = dz_data,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovDZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'X2')),
mxFitFunctionML()),
mxFitFunctionMultigroup(c('MZ','DZ')),
# Output tables
mxAlgebra(A + C + E, name = 'V'),
mxAlgebra(cbind(V, A, C, E),
dimnames = list('X', c('Total', 'A', 'C', 'E')),
name = 'Raw_variance'),
mxAlgebra(cbind(A / V, C / V, E / V),
dimnames = list('X', c('A(%)', 'C(%)', 'E(%)')),
name = 'Variance_components')
)
model <- mxRun(model)
summary(model)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 -0.01540773 0.06578840
## 2 ACE.a[1,1] a 1 1 0.73470965 0.08905756
## 3 ACE.c[1,1] c 1 1 0.56030417 0.12403417
## 4 ACE.e[1,1] e 1 1 0.41025897 0.02890134
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 4 396
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 984.723
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/400
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 192.723 992.723
## BIC: -1113.411 1005.916
## | Sample-Size Adjusted
## AIC: 992.9281
## BIC: 993.2438
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:50:49
## Wall clock time: 0.04244709 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)# Bivariate ACE model
model2 <- mxModel(
name = 'ACE',
# Means
mxMatrix(type = 'Full',
nrow = 1, ncol = 2,
free = TRUE,
name = 'mean'),
mxAlgebra(cbind(mean, mean),
name = 'expMeans'),
# Variance components
mxMatrix(type = 'Lower',
nrow = 2, ncol = 2,
free = TRUE,
name = 'a'),
mxMatrix(type = 'Lower',
nrow = 2, ncol = 2,
free = TRUE,
name = 'c'),
mxMatrix(type = 'Lower',
nrow = 2, ncol = 2,
free = TRUE,
name = 'e'),
mxAlgebra(t(a) %*% a, name = 'A'),
mxAlgebra(t(c) %*% c, name = 'C'),
mxAlgebra(t(e) %*% e, name = 'E'),
# Covariance
mxAlgebra(rbind(cbind(A + C + E, A + C ),
cbind(A + C, A + C + E)),
name = 'expCovMZ'),
mxAlgebra(rbind(cbind(A + C + E, .5%x%A + C),
cbind(.5%x%A + C, A + C + E)),
name = 'expCovDZ'),
# MZ submodel
mxModel(name = 'MZ',
mxData(observed = mz_data2,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovMZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'Y1', 'X2', 'Y2')),
mxFitFunctionML()),
# DZ submodel
mxModel(name = 'DZ',
mxData(observed = dz_data2,
type = 'raw'),
mxExpectationNormal(covariance = 'ACE.expCovDZ',
means = 'ACE.expMeans',
dimnames = c('X1', 'Y1', 'X2', 'Y2')),
mxFitFunctionML()),
mxFitFunctionMultigroup(c('MZ','DZ')),
# Output tables
mxAlgebra(A + C + E, name = 'V'),
mxAlgebra(abs(A) + abs(C) + abs(E), name = 'Vabs'),
mxAlgebra(abs(A) / Vabs, name = 'Aperc'),
mxAlgebra(abs(C) / Vabs, name = 'Cperc'),
mxAlgebra(abs(E) / Vabs, name = 'Eperc'),
mxAlgebra(cbind(diag2vec(Aperc),
diag2vec(Cperc),
diag2vec(Eperc)),
dimnames = list(c('X', 'Y'), c('A(%)', 'C(%)', 'E(%)')),
name = 'Variance_components'),
mxAlgebra(cbind(vech(V),
vech(A), vech(C), vech(E),
vech(Aperc), vech(Cperc), vech(Eperc)),
dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
c('Total', 'A', 'C', 'E', 'A(%)', 'C(%)', 'E(%)')),
name = 'Covariation_components'),
mxMatrix(type = "Iden",
nrow = 2, ncol = 2,
name = "I"),
mxAlgebra(sqrt(I * V), name = "SD_total"),
mxAlgebra(sqrt(I * A), name = 'SD_A'),
mxAlgebra(sqrt(I * C), name = 'SD_C'),
mxAlgebra(sqrt(I * E), name = 'SD_E'),
mxAlgebra(solve(SD_total) %&% V, name='r_total'),
mxAlgebra(solve(SD_A) %&% A, name = 'r_A'),
mxAlgebra(solve(SD_C) %&% C, name = 'r_C'),
mxAlgebra(solve(SD_E) %&% E, name = 'r_E'),
mxAlgebra(cbind(vech(r_A), vech(r_C), vech(r_E),
vech(r_total)),
dimnames = list(c('X <-> X', 'X <-> Y', 'Y <-> Y'),
c('r_A', 'r_C', 'r_E', 'Total')),
name='All_correlations')
)
model2 <- mxRun(model2)
summary(model2)
## Summary of ACE
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 ACE.mean[1,1] mean 1 1 0.028582935 0.06602392
## 2 ACE.mean[1,2] mean 1 2 -0.006697959 0.06398142
## 3 ACE.a[1,1] a 1 1 0.229626300 0.19038682
## 4 ACE.a[2,1] a 2 1 0.487155814 0.14554958
## 5 ACE.a[2,2] a 2 2 0.761981284 0.09837065
## 6 ACE.c[1,1] c 1 1 0.525659559 0.07115918
## 7 ACE.c[2,1] c 2 1 0.367418707 0.20941586
## 8 ACE.c[2,2] c 2 2 0.456054195 0.15874248
## 9 ACE.e[1,1] e 1 1 0.417843143 0.02948765
## 10 ACE.e[2,1] e 2 1 0.493956895 0.05375951
## 11 ACE.e[2,2] e 2 2 0.467266025 0.03298568
##
## Model Statistics:
## | Parameters | Degrees of Freedom
## Model: 11 789
## Saturated: NA NA
## Independence: NA NA
## | Fit (-2lnL units)
## Model: 1793.504
## Saturated: NA
## Independence: NA
## Number of observations/statistics: 200/800
##
## Information Criteria:
## | df Penalty | Parameters Penalty
## AIC: 215.5039 1815.504
## BIC: -2386.8685 1851.785
## | Sample-Size Adjusted
## AIC: 1816.908
## BIC: 1816.936
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2021-01-12 03:50:49
## Wall clock time: 0.06305766 secs
## optimizer: SLSQP
## OpenMx version number: 2.18.1
## Need help? See help(mxSummary)mxEval(Covariation_components, model2)
## Total A C E A(%)
## X <-> X 1.1199498 0.2900490 0.4113145 0.4185863 0.2589839
## X <-> Y 0.7695757 0.3712036 0.1675628 0.2308093 0.4823484
## Y <-> Y 1.0069384 0.5806155 0.2079854 0.2183375 0.5766147
## C(%) E(%)
## X <-> X 0.3672615 0.3737545
## X <-> Y 0.2177341 0.2999176
## Y <-> Y 0.2065523 0.2168331library(mlth.data.frame)
# mlth.data.frame - multiheaded data.frame
model2 <- multivariate_ace(
data = data2,
zyg = 'zyg',
ph = c('X', 'Y')
)
model2 <- def_ci(model2, ci = 'Variance_components')
model2 <- mxRun(model2, intervals = TRUE)
get_output_tables(model2)
## $Variance_components
## A(%)-------------------------- C(%)---------------------------- E(%)-------------------------
## est lCI uCI est lCI uCI est lCI uCI
## X 0.2589905 0.02667819 0.6127657 0.3672559 4.764229e-02 0.5861756 0.3737536 0.2791622 0.4919170
## Y 0.5766131 0.31274212 0.8272015 0.2065536 5.511957e-33 0.4481345 0.2168333 0.1594465 0.2973594
##
## $Covariation_components
## Total A C E A(%)
## X <-> X 1.1199492 0.2900561 0.4113080 0.4185850 0.2589905
## X <-> Y 0.7695746 0.3712059 0.1675600 0.2308088 0.4823520
## Y <-> Y 1.0069377 0.5806135 0.2079866 0.2183376 0.5766131
## C(%) E(%)
## X <-> X 0.3672559 0.3737536
## X <-> Y 0.2177306 0.2999174
## Y <-> Y 0.2065536 0.2168333
##
## $All_correlations
## r_A r_C r_E Total
## X <-> X 1.0000000 1.0000000 1.0000000 1.0000000
## X <-> Y 0.9045451 0.5728871 0.7634776 0.7246867
## Y <-> Y 1.0000000 1.0000000 1.0000000 1.0000000# Formatting the table for html or latex report
# Powered by knitr and kableExtra
library(kableExtra)
M %>%
behead() %>%
kable(digits = 3,
align = 'c') %>%
add_complex_header_above(M)|
A(%)
|
C(%)
|
E(%)
|
|||||||
|---|---|---|---|---|---|---|---|---|---|
| est | lCI | uCI | est | lCI | uCI | est | lCI | uCI | |
| X | 0.259 | 0.027 | 0.613 | 0.367 | 0.048 | 0.586 | 0.374 | 0.279 | 0.492 |
| Y | 0.577 | 0.313 | 0.827 | 0.207 | 0.000 | 0.448 | 0.217 | 0.159 | 0.297 |
# Saving the table for the output
M %>%
register_output(
name = 'Table 1',
caption = 'Table 1: Variance components')
## A(%)-------------------------- C(%)---------------------------- E(%)-------------------------
## est lCI uCI est lCI uCI est lCI uCI
## X 0.2589905 0.02667819 0.6127657 0.3672559 4.764229e-02 0.5861756 0.3737536 0.2791622 0.4919170
## Y 0.5766131 0.31274212 0.8272015 0.2065536 5.511957e-33 0.4481345 0.2168333 0.1594465 0.2973594